/*
 * EllipsePore.cpp
 *
 *  Created on: Nov 14, 2011
 *      Author: yi
 */

#include "mbuild/EllipsePore.h"

EllipsePore::EllipsePore() {
	myfunc=NULL;
	funcdensitydist=NULL;
	holefunc=NULL;
}

EllipsePore::~EllipsePore() {
}

EllipsePore::EllipsePore(double R_in, double R_out, double T)
{
	setup(R_in,R_out,T);
}

void EllipsePore::setup(double R_in, double R_out, double T)
{
	this->R_in=R_in;
	this->R_out=R_out;
	if(R_in >= R_out) throw ParameterError();
	this->T=T;
	if(T <= 0) throw ParameterError();
}



inline double EllipsePore::f(double x, double y, double z)
{
	double r=sqrt(x*x+y*y);
	if(r<R_out) return (r-R_out)*(r-R_out)/((R_out-R_in)*(R_out-R_in))+4*z*z/(T*T)-1;
	else if (z > 0) return z-T/2;
	else return -T/2-z;
}



double EllipsePore::zf(double x, double y)
{
	double r=sqrt(x*x+y*y);
	double z=1-(r-R_out)*(r-R_out)/((R_out-R_in)*(R_out-R_in));
	z=sqrt(T*T*z/4);
	if (r<= R_out) return z;
	else return T/2;
}



inline double EllipsePore::df_dx(double x, double y, double z)
{
	double r=sqrt(x*x+y*y);
	if(r<=R_out) return 2*(r-R_out)*x/(r*(R_out-R_in)*(R_out-R_in));
	else return 0;
}



inline double EllipsePore::df_dy(double x, double y, double z)
{
	double r=sqrt(x*x+y*y);
	if(r<=R_out) return 2*(r-R_out)*y/(r*(R_out-R_in)*(R_out-R_in));
	else return 0;
}





inline double EllipsePore::df_dz(double x, double y, double z)
{
	double r=sqrt(x*x+y*y);
	if(r<=R_out) return 8*z/(T*T);
	else if (z >0 )return 1;
	else return -1;
}



